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<N . Abstract 

a 

We include the effect of neutrino free streaming into the spectrum of relic gravitational waves (RGWs) 
^ . in the currently accelerating universe. For the realistic case of a varying fractional neutrino energy 

density and a non-vanishing derivative of mode function at the neutrino decoupling, the integro- 



CO 



differential equation of RGWs is solved by a perturbation method for the period from the neutrino 
decoupling to the matter-dominant stage. Incorporating it to the analytic solution of the whole history 
of expansion of the universe, the analytic solution of GRWs is obtained, evolving from the inflation 



> 
O 
m t 

up to the current acceleration. The resulting spectrum of GRWs covers the whole range of frequency 
(10~ 19 ~ 10 10 )Hz, and improves the previous results. It is found that the neutrino free-streaming 
causes a reduction of the spectral amplitude by ~ 20% in the range (10~ 16 ~ 10~ 10 ) Hz, and leaves the 
other portion of the spectrum almost unchanged. This agrees with the earlier numerical calculations. 
Examination is made on the difference between the accelerating and non-accelerating models, and our 



analysis shows that the ratio of the spectral amplitude in accelerating ACDM model over that in CDM 
model is ~ 0.7, and within the various accelerating models of > Q, m the spectral amplitude is 



proportional to $7 m /SlA for the whole range of frequency. Comparison with LIGO S5 Runs Sensitivity 
shows that RGWs are not yet detectable by the present LIGO, and in the future LISA may be able to 
detect RGWs in some inflationary models. 

PACS numbers: 04.30.-w, 98.80.-k, 04.62,+v 

1. Introduction 

Inflationary models generally predict the existence of a stochastic background of relic gravitational waves 
(RGWs). [1, 2, 3]. Due to their very weak coupling with matter, RGWs still encode a wealth of information 
about the very early universe when they were generated, and enable us to study the inflationary and the successive 
physical processes, much earlier than the recombination time at a temperature T ~ 0.3 ev, up to which CMB 
information can tell. Not only the RGWs are the scientific goal of the detections, such as the laser interferometers 
now underway [4, 5, 6, 7], but also are a source, along with the density perturbations, of CMB anisotropics and 
polarizations [8, 9, 10, 11, 12, 13]. In particular, the B-polarization of CMB can only be generated by RGWs. 
Thus, it is important to calculate the spectrum of RGWs, which depends on several physical processes. First of all, 
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RGWs can be further modified by the subsequent expansion of the universe, giving rise to the redshift-suppression 
on the spectrum. In our previous analytic and numerical investigations [3], we studied the RGWs in the current 
accelerating expansion of the universe, obtained the modifications on the spectrum by the presence of dark energy. 
In particular, we have found that the amplitude of RGWs is reduced by a factor ~ 0.3 in comparison with the 
matter-dominant models, and that within the ACDM models with > ^m, the amplitude oc Q, m /Q^, over 
almost the whole frequency range of the spectrum. 

There are other processes that can also change the spectrum of RGWs. One important process is the free- 
streaming of neutrinos that occurred in the early universe [14]. It will leave the imprints on the spectrum. At a 
temperature T ~ 2 Mev during the radiation-dominant stage in the early universe, cosmic neutrinos decoupled 
from electrons and photons, and started free-streaming in space. This will give rise to an anisotropic part 7Ty 
of the energy-momentum tensor Tij as a source of the equation of RGWs, and will cause a damping effect on 
the RGWS. Weinberg analyzed the effect and arrived at the integro-differential equation for RGWs, and gave an 
estimate of the damping on RGWs due to the neutrinos free-streaming [14]. Subsequently, in the special case 
of a constant fractional neutrino energy density /^(0), a vanishing time u^ ec = of the neutrino decoupling, 
and a vanishing time-derivative x'{ u dec) = 0, Dicus and Repko [15] obtained an analytic solution, in terms of a 
series of Bessel's functions, of the integro-differential equation for the radiation stage, qualitatively agreeing with 
Weinberg's estimate. However, this solution holds only for the short wavelength modes reentering the horizon 
long after the neutrino decoupling during radiation-dominant stage, and the conditions it has used are actually 
approximations and will obviously cause some errors. Moreover, as Weinberg points out, the solution for the 
radiation stage is still to be joined with those for other expansion stages, so that the effect of neutrino free 
streaming is taken into account in a complete computation of the spectrum of RGWs. In a numerical study of the 
matter-dominant universe, Watanabe and Komatsu [16] investigated the damping effects on the RGWs caused 
by the evolution of the effective relativistic degrees of freedom, including the neutrino free-streaming, and gave a 
numerical solution of the energy density spectrum [16]. But there the important effect of the acceleration of the 
present universe has not been considered. 

In this paper, extending our previous work on the analytic spectrum of RGWs in the accelerating universe 
[3], we will include the damping effect of neutrino free-streaming into our analytic calculation scheme. Both 
effects of the accelerating universe and of the neutrino free streaming are taken into account, simultaneously. The 
following improvements are achieved in this paper over the previous studies. In comparison with Ref.[16], the 
damping effect on the spectrum of RGWs by the dark energy Q\ is now properly included. Different from Dicus 
and Repko's method of series expansion that is valid in the special case [15], we apply a perturbation method 
to solve the integro-differential equation of RGWs by an iterative procedure. Actually, for practical use, the first 
order solution is enough for an evaluation for the spectrum of RGWs, and the solution of higher accuracy can 
be easily achieved by going to higher order of iterations. This calculation has the merit of precisely taking into 
account of the time-varying fractional neutrino energy density f u {u), the non-vanishing time u& ec / and the 
non-vanishing time-derivative x'( u dec) 7^ of the mode functions. Therefore, the result is valid for all the modes 
of RGWs of an arbitrary wavelength, and reduces to that in Ref. [15] in the special case for the short wave length 
limit. We give the analytic expressions of the full spectrum h(k,r]H) of RGWs itself and of the spectral energy 
density Q g (k), valid for the whole range of frequencies. As a comprehensive compilation, by using the parameters 



(3, (3 S , 7 and r, respectively, such important cosmological elements have been explicitly parameterized, as the 
inflation, the reheating, the dark energy, the tensor/scalar ratio. This will considerably facilitate further studies on 
the RGWs and the relevant physical processes. Besides, several typographical errors in the previous studies have 
been corrected thereby. So not only can it be easily used in computation for other applications in cosmology, such 
as calculations of CMB anisotropics and polarizations generated by RGWs [12], but also can be directly compared 
with the sensitivity curves of those ongoing and forthcoming laser interferometer GW detectors, such as LIGO, 
LISA, etc [4, 5, 6, 7]. 

The outline of this paper is as follows. In section 2, to various stages of expansion of the universe the scale 
factor a(rj) is specified with the parameters being determined by the continuity conditions. In section 3, we present 
the analytical solutions of the modes hk(rj) of RGWs during each stage, and, in particular, include the effect of 
neutrino free streaming during the radiation-dominant stage. The subtleties of interpreting the observational data 
within the non-accelerating models are discussed. In section 4, we present the resulting spectrum and analyze the 
effects of (3, (3 S , 7, r and the neutrino free streaming. The Appendix gives the detailed calculation of the anisotropic 
part 717, of the energy-momentum tensor, and present the perturbation method for the solution modified by the 
neutrino free-streaming during the radiation stage. In this paper we use unit with c = U = ks = 1. 

2. Expansion history of the universe 

From the inflationary up to the current accelerating stage, the expansion of a spatially flat universe can be 
described by the spatially flat (Q\ + £l m + Q r = 1) Robertson-Walker spacetime with a metric 

ds 2 = a 2 (r])[-dr] 2 + SijdjdaP], (1) 

where the scale factor has the following forms for the successive stages [17]: 
The inflationary stage: 

a(v) = lo\v\ 1+(S , -oo< V < m , (2) 

where 1 + (3 < 0, and 771 < 0. The special case of (3 = — 2 is the de Sitter expansion of inflation. 
The reheating stage: 

a{v) = a z \v ~ Vp\ 1+P % m<V<Vs, (3) 

here we take the absolute value of rj — rj p , different from Ref. [3]. This is because 1 + (3 S might be negative for 
some models of the reheating. As a model parameter, we will mostly take (3 S = —0.3, though other values are 
also taken to demonstrate the effect of the various reheating models. 
The radiation-dominant stage: 

a(rj) = a e (rj - rj e ), r/ s < rj < 7/2. (4) 

This is the stage during which the neutrinos decoupled from the radiation component. We use r/rf ec to denote 
the starting time of the neutrino decoupling: r] s < r]d ec < i]2- The corresponding energy scale is ~ 2 Mev for the 
decoupling. As will be seen later, the wave equation of RGWs is still homogeneous for r\ < r]d ec , but becomes 
in homogeneous for r/^ec < V < V2- 
The matter-dominant stage: 

a(n) = a m (v ~ Vm) 2 , V2<V<VE- (5) 



The accelerating stage up to the present time rfn [3]: 

a{v) = Ih\v ~ %r 7 , Ve < V < W, (6) 

where the index 7 depends on the dark energy Q A . By numerically solving the Friedmann equation [3], 

a' 8ttG 

(-2) = -5- (PA + Pm + p r ), (7) 



where a' = da/dr/, we find that 7 ~ 1.06 for n A = 0.65, 7 ~ 1.05 for Q A = 0.7, and 7 ~ 1.044 for tt A = 0.75 
(as a correction to 7 ~ 1.048 in Ref.[3]). 

In the above specifications of a(rj), there are five instances of time, 771, rj s , 772, tje, and T]h, which separate the 
different stages. Four of them are determined by how much a(rj) increases over each stage by the cosmological 
considerations. We take the following specifications: Ci = = 300 for the reheating stage, ( s = ^np; _ iq 24 
for the radiation stage, C2 = ^P4 = ^4^\ = 3454C 1 for the matter stage, and Ce = 4^ = (^) 1/3 for 
the present accelerating stage. Note that here i^) 1 ^ is model-dependent, and associated with the value of 7, 
instead of the fixed value (1.33, as in Ref.[3]). The remaining time instance is fixed by an overall normalization, 
namely 

\VH ~ Va\ = 1- (8) 

There are twelve constants in the expressions of a(rj), among which (3, (3 S and 7 are imposed as the model 
parameters, for the inflation, the reheating, and the acceleration, respectively. So there remain nine constants. By 
the continuity of a(rj) and of a(rf)' at the four given joining points 771, 7? s , 772 and 77^, one can fix eight constants. 
Only one constant remains, which can be fixed by the present expansion rate Ho of the universe and Eq.(8), 

l H = 7/H . (9) 

Then all parameters are fixed as the following: 

2 i 

VE~Vm = ~Q, 

7 

2 — i - 

V2-Vm = -C 2 2 Q, 

7 

1 -I A 
V2~Ve = -C 2 2 Q, 

7 

1 -i - 

Vs~Ve = -C S _1 C 2 2 Ce> 

7 

77 s -77 P = i(l + /3 s )C 1 C 2 _l ci 
7 

m - v P = ^(1 + /? s )Ci T ^ 7 C 1 C 2 ~"cl, 

7/ dec = 1.15 x 10- 10 C£C 2 r? 2 , (10) 



and 

a m = -7- 7^ Ck 



a e = Ih 7 C 2 Ce > 

a, = Zh7 1+/J '|1 + A|- (1+A) C^cf^C^ Cl+i ^ ) , 

io = iH7 1+ "\i + Pr {1+ " ) cr s ck^c E T } . (ii) 

The above expressions correct some typographical errors in Ref.[3]. 

In the expanding universe, the physical wavelength is related to the comoving wavenumber k by 

A^^M (12) 

and the wavenumber kn corresponding to the present Hubble radius is 

27ra(r? H ) 

= i/ Hq = ( 13 ) 

There is another wavenumber 

_2tto(j]e) k H 

kE = ~T/h^ = Tw (14) 

whose corresponding wavelength at the time tje is the Hubble radius 1/Hq. In the present universe the physical 
frequency corresponding to a wavenumber k is given by 

A 2ira(r] H ) 27r 7 v ' 

3. Analytical solution 

In the presence of the gravitational waves, the perturbed metric is 

ds 2 = a 2 (r/)[-dr/ 2 + (<% + h^dx^x^, (16) 
where the tensorial perturbation ^ is a 3 x 3 matrix and is taken to be transverse and traceless 

h\ = 0, hijj = 0. (17) 

The wave equation of RGWs is 

d„(V=gd»h tj ) = 0. (18) 

However, from the temperature T ~ 2 Mev up to the beginning of the matter domination, the neutrinos are 
decoupled from electrons and photons and start to freely stream in space. This effect of neutrino free streaming 
gives rise to an anisotropic portion irij of the energy-momentum stress Tij. Then Eq.(18) acquires an inhomo- 
geneous source term —16TrGirij on the right hand side during the period rjdec < V < As is shown in the 
Appendix, the anisotropic stress is also transverse and traceless, and it is zero before the decoupling and 
becomes negligible small after the matter domination. To solve the equation, we decompose hij into the Fourier 
modes of the comoving wave number k and into the polarization state a as 

d 3 k 



Mifcx) = E/ 7j^44%K kx , (19) 



where h^l*(r)) = h^\rf) ensuring that be real, e?- is the polarization tensor, and a denotes the polarization 
states x,+. Here hij is treated as a classical field, instead of a quantum operator [1, 3]. In terms of the mode 
h ( k a \ Eq.(19) reduces to 

4 CT %) + 2^4 CT)/ W + ^ 2 4 CT) W = 0. (20) 



Since for each polarization, x, +, the wave equation is the same and has the same statistical properties, from 
now on the super index (a) can be dropped from h k a \ As demonstrated in Eq.(2) through Eq.(6), for all the 
stages of expansion the time-dependent scale factor is of a generic form 

a(v) oc r] a , (21) 
the solution to Eq.(20) is a linear combination of Bessel function J u and Neumann function N v 

h k ( V ) = x l 2- a [ ai J Q _i (k V ) + a 2 N a _i(k V )}, (22) 

where the constants a\ and a 2 are determined by the continuity of h k and of h' k at the joining points r]i,r) s ,r) 2 
and t]e- However, as mentioned earlier, during the neutrino free streaming with r)d ec < f] < i]2, Eq.(20) will be 
modified and its solution will be given later. 
The inflationary stage has the solution 

h k {r]) = A Q l l \r ] \-^-t 3 [A l Ji_ +f3 {x) + A 2 J_ ( i +p) {x)), -(xXr]< m (23) 

where x = krj and 

A x = 1 A 2 = iA ie -^, (24) 

COS (jTT V 2 

are taken [18], so that the so-called adiabatic vacuum is achieved: lim^oo h k (r/) oc e~ tkv in the high frequency 
limit [19]. Moreover, the constant Aq in Eq.(23) is independent of k, whose value is determined by the initial 
amplitude of the spectrum, so that for kr] <C 1 the fc-dependence of hk(i]) is given by 

h k (v) oc Ji+p(x) « (25) 

As will be seen, this choice will lead to the required initial spectrum in Eq.(45). 
The reheating stage has 



h k (v) = t~~ Ps BiJi+^kt) + B 2 N l+/3s (kt) , m<V<Vs (26) 



where t = rj — i] p and 



B 1 = -^irtf^'ikN^iktJhkirn) + N^k^h'M], (27) 
B 2 = ^4 +Ps [kJ l+Ps (kt 1 )h k {r h ) + Ji+^ktMm)}, (28) 

with ti = rji — i] p , and h k (r]i) and h' k (r]i) are the corresponding values from the precedent inflation stage. 

The radiation-dominant stage needs to be divided into two parts. The first part of the stage is before the 
neutrino decoupling when r] s < i] < r]d ec , the neutrino damping is ineffective yet, the wave equation is still 
homogenous with the solution 



1 



where y = i] — rj e and 



h k (n) = y ^[CiJi(ky) + C 2 Ni(ky)\, r, s < rj < r] dec (29) 



Ci = -^iryl[kN3(ky s )h k (ri s ) + Ni(ky 8 )h' k (rj 9 )], (30) 
1 ± 

C 2 = -nyl[kJz(ky s )h k (r) s ) + Ji(ky s )ti k (r] s )], (31) 



where y s = rj s — rf e , and hk{r} s ) and h' k (rj s ) are from the reheating stage. If we do not include the neutrino effect 
and let Eq.(29) be valid for the whole radiation stage rf s < 77 < 7/2, then our previous exact result [3] would be 
recovered . 

The second part is from the neutrino decoupling up to the matter domination with r\dec < V < V2- The 
temperature at the neutrino decoupling r/^ec is taken to be T ~ 2 Mev. During this period the wave equation is 

K(V) + 2 ^ h 'k(n) + k 2 h k ( V ) = levrGa 2 ^^). (32) 
The detailed derivation are given in the Appendix [14]. Writing the mode function as 

h k (v) = h k (Vdec)x(v), Vdec <V<V2 (33) 

where /ifc(%ec) is given by Eq.(29) evaluated at r^ec. and x( u ) satisfies the following integro-differential equation 

X» + -X'(u) + X (u) = -24 f{ 0) f dUK(u - U) X '(U), (34) 
u u\l + au) J Udec 

where u = krf, /^(O) = 0.40523 is the fractional energy density of neutrinos at u = 0, a = a e /(k 0(772)) = 
fc^^C^Ci/ 7 . and K(u) is the kernel defined in Eq.(72) in the Appendix. In dealing with Eq.(34) Dicus & Repko 
[15] use the following approximations 

a = 0, u dec = 0, x'iVdec) = 0, (35) 

and derive an analytical solution, valid for those modes reentering the horizon long after the neutrino decoupling 
during the radiation-dominant stage. Here without making the approximations in Eq.(35), we try to give a solution 
valid for all the modes that reenter the horizon both before and after the neutrino decoupling. The idea is that, 
the neutrino damping effect is small, therefore the right hand side of Eq.(34) will cause only a small variation 
to the homogeneous solution Xo(^)- Thus, as an approximation, one substitutes Xo{ u ) i n place of x( u ) in t ne 
integration on the right hand side of Eq.(34), and obtains the first order approximate solution. As it turns out, 
this is accurate enough for our purpose of calculating the spectrum. For the higher order solutions this process can 
be iterated to achieve higher accuracy. The Appendix gives the detailed expressions of the solutions. To examine 
this perturbation method, Figure 1 plots the mode function x( u ) as the solution of different orders, respectively, 
where the same assumption as Eq.(35) are adopted to compare with Dicus ii Repko's analytic result [15]. It is 
seen that the first, and second order solutions differ by ~ 4%, and by ~ 1%, respectively, from the exact one, the 
third order solution almost overlaps with it. Therefore, our method is effective in evaluating the damping caused 
by neutrino free-streaming. Moreover, our result, Eqs.(83) and (86) in the Appendix, holds for any wavelength 
and for the realistic condition of a ^ 0, u& ec / 0, and x'iVdec) / 0, while the approximation in Refs.[14] and 
[15] is not valid for the very short nor for the long modes. Figure 1 shows that the solution without the neutrino 
free-streaming has a higher amplitude, as expected. 

Our calculation reveals that the neutrino damping on the RGWs is mainly pronouncing only in the frequency 
range v ~ (10~ 16 ~ 10~ 10 ) Hz, which corresponds to k ~ (10 2 ~ 10 8 ) and a ~ (10~ 7 ~ 10" 1 ). Outside this 
range the neutrino damping barely alters the RGWs. Figure 2 shows that our solution of short (u > 10~ 10 Hz) 
and long {y < 10~ 16 Hz) modes almost overlap the homogeneous solution without neutrino free-streaming. For 

thp <;hnr1- mnHp? rppntprincr thp hnri7nn wpII hpfnrp thp Hprnnnlincr thp fartnr 1 lip" nn thp r h <; nf Fn (7R\ is 
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Figure 1: Under the approximations in Eq.(35), the solutions by our method are compared with that in 
Ref.[15]. 




Figure 2: Neutrino free-streaming barely affects the short (upper) and the long modes (lower). 

very small and the inhomogeneous term is negligible. For those long modes, they are still outside the horizon 
during the neutrino free streaming, and are not affected by the damping. Only much later do these modes 
reenter the horizon, the neutrino density f v {u) becomes negligibly small, the homogeneous solution is valid for 
these long modes. Therefore, in long and short wavelength limit, the solution for RGWs is practically that of 
the homogeneous equation. Moreover, the upper panel of Figure 2 shows that, at the neutrino decoupling time 
rjdec, the time derivative of the short mode function x'iVdec) deviates from zero considerably. Therefore, the 
approximation in Eq.(35) is not accurate enough for the short modes of RGWs. 
The matter-dominant stage has 



3 



h k (r])=z-i D 1 J 3 (kz)+D 2 N 3 (kz) , m < V < Ve (36) 

L 2 2 J 



where z = rj — rj m and 



Di = —it zl[kNs{k z 2 )h k { m ) + Ns{k z 2 )h' k { m )l (37) 

1 £ 

D 2 = -vrz| [kJs,{kz 2 )h k (r] 2 ) + Jz (k z 2 )h' k (r] 2 )], (38) 

2 2 2 

withe z 2 = 7] 2 — r\m- In the expressions of D\ and D 2 , the mode functions h k (r] 2 ) and h' k (r] 2 ) are again from the 
nrpredpnt ^tacrp 
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Figure 3: The scale factor a(t) in the accelerating, and non-accelerating models, respectively. 

The accelerating stage has 

h k (v) =s^[E 1 J_ 1 __^ks) + E 2 N_ 1 _^(ks)\, V E<V<VH (39) 

where s = 77 — rj a and 

El = ^sl^lkN^iksE^M+N^iksE^ivE)], (40) 
E 2 = ^s E ~ J [kJi^{ks E )h k {rj E ) + J_i_^(k8 B )ti k (r]E)], (41) 

with se = f]E — Va- So far, the explicit solution of hk{rj) has been obtained for all the expansion stages, from 
Eq.(23) through Eq.(39). 

The above detailed expressions of h k {rj) are the major ingredients to determine the the spectrum of RGWs 
in the accelerating universe. What kind of RGWS would a matter-dominant universe have? To compare with 
the spatially flat accelerating universe, this non-accelerating universe is assumed to be also spatially flat with 
Q m + O r = 1. It should also go through the consecutive expansion stages listed previously, from the inflation 
to the matter-dominant, except the accelerating stage that is replaced by a continuation of the matter-dominant 
stage up to the present time tjh- In each stage the mode function hk(rj) is of the same form as those given in 
Eqs.(23), (26), (29), (33), and (36), respectively. But the time duration of the matter stage for Eq.(36) is now 
extended to 7/2 < rj < t]h- In both the accelerating and the matter-dominant models the mode hk(rf) is sensitive 
to the scale factor a(rj) determined by their respective Friedmann equation (7), in which one sets p\ = for 
the matter-dominant model. To have a specific comparison of the two models, let us start at the time 772 of the 
equality of radiation-matter with 1 + z = 3454, when pa « /> ra = p r , and it can be assumed that both models 
have the same initial values 0(772) and a' (772). Instead of the sudden transition approximation as in Eqs.(5) and 
(6), we solve numerically the Friedmann equation in both models up to the present time tjh with z = 0. Doing this 
is equivalent to assuming that both models would have an equal age of the universe. As a result, it is found that 
the scale factor a(rjH) i n the accelerating model is ~ 1.3 times of that in the matter-dominant model shown in 
Fig. 3. This difference of a(r)H) will consequently cause a difference in the spectra for the two models. As is known 
fl. 31. inside the horizon the arrmlitude of modes hu(n) oc l/a(n). so the matter-dominant model would Dredict 



a spectral amplitude higher than the accelerating model. Indeed, our analytic calculation demonstrates that the 
ratio of the spectral amplitudes of CDM over those of ACDM is ~ 1.3. Moreover, we like to emphasize that 
there are some subtleties with the matter-dominant models, regarding to interpretation of the current cosmological 
observations. The actual universe is an accelerating one, so the observed Hubble constant is properly interpreted 
as the current expansion rate in the accelerating model, Hq = a' /a 2 (r]H)- However, as our calculation has shown, 
the virtual matter-dominant universe would have a smaller rate a' /a 2 (r]H) — 0.65 Hq. In this regard, Ref.[16] uses 
the observed Hubble constant Hq as the current expansion rate of the virtual matter-dominant universe. This 
would give a spectrum with amplitude lower by an extra factor ~ 1.3 than it should have. 

4. Spectrum of relic gravitational waves 

The spectrum of RGWs h(k,rf) at a time rj is defined by the following equation [17]: 



re 

J /ri/,-.//;^ U/,^x./ / )/,;,-lx.//)() . (12) 



where the right-hand side is the expectation value of the h^hij. Calculation yields the spectrum, which is related 
to the mode function h k {rj) as follows 

h(k, V ) = -k^ 2 \h k (ri)\, (43) 

7T 

where the factor 2 counts for the two independent polarizations. At present with time r\n the spectrum is 

h(k, m ) = -k 3 / 2 \h k ( m )\. (44) 

7T 

Note that this expression is formally different from the previous one in Refs.[3] only because here we use a different 
expansion for hij(r), x) in Eq.(19). One of the most important properties of the inflation is that the initial spectrum 
of GRWs at the time rji of the horizon-crossing during the inflation is nearly scale-invariant [17]: 

h(k,T H ) = A(-£-) 2 +P, (45) 

where 2 + (3 ~ 0, and A is a /c-independent constant to be fixed by the observed CMB anisotropics in practice. 
The First Year WMAP gives the scalar spectral index n s = 0.99 ± 0.04 [9]. The Three Year WMAP gives 
n s = 0.95ljlooi9 [10] > while in combination with constraints from SDSS, SNIa, and the galaxy clustering, it 
would give n s = 0.965 ± 0.012 (68% CL) [11]. From the relation n s = 2(3 + 5 [1, 3], we have the inflation 
index = —2.02 for n s = 0.951. Note that the constant A is directly proportional to Aq in Eq.(23) through the 
relation (43). Since the observed CMB anisotropies [9] is AT/T ~ 0.37 x 10~ 5 at / ~ 2, which corresponds to 
anisotropies on scales of the Hubble radius 1/Hq, so, as in Refs.[3], we take the normalization of the spectrum 

h(k E ,7] H ) = 0.37 x W~ 5 A , (46) 

where k E = (i+" E ) = (i+z E ) ' s t ^ ie wave number that crosses the horizon at tje, its corresponding physical 
frequency being ve = kE/2ira(r]H) = Hq/(1 + z e ) ~ 1 x 10~ 18 Hz, r is taken as a parameter roughly representing 
the tensor/scalar ratio. In Eq.(46) it is r5 rather than r as in Ref.[3]. The value of the ratio r is an important 
issue and is still unsettled yet. However, as examined in details in Ref.[20], the relative contributions from the 
RGWs and from the density perturbations are, in fact, frequency-dependent; thus, generally speaking, for different 



frequency ranges r can take on a different values. Therefore, in our treatment, for simplicity, r is only taken as a 
constant parameter for normalization of RGWs, and does not accurately represent the actual relative contributions. 
Currently, only observational constraints on r have been given. Recently the Three Year WMAP constraint is 
r < 2.2 (95% CL) evaluated at k = 0.002 Mpc" 1 , and the full WMAP constraint is r < 0.55 (95% CL) [12]. 
The combination from such observations, as of the Lyman-a forest power spectrum from SDSS, 3-year WMAP, 
supernovae SN, and galaxy clustering, gives an upper limit r < 0.22 (95% CL) [11]. Moreover, the ratio r may 
be allowed to take on different values on different range of frequency, but we will take a constant r for simplicity. 
The spectral energy density tt g (k) of the RGWs is given by 

2 7 

n g (k) = ^-h 2 (k, VH )(—) 2 , (47) 

directly associated with the spectrum of RGWs h(k,rjH) in Eq.(44). This follows from the definition [17, 20] 

Pg fkupper dk 

Mgw = — = to g (k) — , (48) 

Pc J how k 

where p g = 32^5^.7,0^0 ' s tne energy density of RGWs, and p c = 3Hq/8itG is the critical energy density. The 
integration in Eq.(48) has the lower and upper limits, ki ow and k upper , as the cutoffs of the wavenumber. For the 
lower limit ki ow , the corresponding wavelength may be taken to be the current Hubble radius, A/ ow = 1/Hq. This 
is because the waves with wavelengths longer than 1/Hq should be treated as part of the space-time background 
and should not be included to the energy of RGWS [1] [21]. By Eq.(15), the corresponding frequency 

v low ~ 2 x 10- 18 Hz. (49) 

The upper limit k upper can be determined by as the following. During the inflation the modes of GRWs with 
wavenumbers greater than the expansion rate H(rji) are approaching the adiabatic limit, therefore, their generation 
is thus effectively suppressed [19]. Taking the scale of the vacuum energy driving the inflation to be E vac ~ 10 16 
Gev, typical of Grand Unified Theories, then H(r)i) ~ 10 13 Gev ~ 10 38 Hz. During the subsequent stages of 
cosmic expansion, the corresponding frequency v of this value will be redshifted by a factor a(rji) / a(rjH) ~ 10~ 29 ; 
thus one has 

"upper ^ 10 10 Hz. (50) 

If the energy scale for the inflation is lower than 10 16 Gev, then the upper limit v upper will be lower than that in 
Eq.(50) correspondingly. These lower and upper integration limits in Eqs.(49) and (50) also ensure the convergence 
of the integration of Eq.(48). 

In the absence of direct detection of RGWs, the constraints on the energy density £lcw is more relevant. 
Given the model parameters (3, (3 S , 7, and r, the definite integration of Eq.(48) yields Q-gw of RGWs. For the 
fixed parameters r = 0.22, Q\ = 0.75, and (3 S = —0.3, one finds Qgw = 1-12 x 10 -2 for the inflationary model 
of (3 = —1.8. Such a large energy density will inevitably affect the expansion rate of the universe at a temperature 
T ~ a few MeV when the nucleosynthesis process is going on. The nucleosythesis bound is [22] 

n GW h 2 < 8.9 x 10~ 6 (51) 

with h ~ 0.71 being the Hubble parameter [9]. Thus the (3 = —1.8 model with r = 0.22 predicts an energy density 
heintr mmp fnur orders hicrher than the unner hnund trivpn Fn TR1 ^ Sn thiq model will he in iennardv nnle^ 




Figure 4: The spectrum h(i/,r)H) of GRW is very sensitive to the inflation parameter (3. 



the parameter r is much smaller than 0.22. Under the same set of parameters r = 0.22, Q\ = 0.75, and j3 8 = —0.3, 
the model of (3 = -1.9 gives Q GW = 2.04 x 10~ 8 , and the model of (3 = -2.02 gives Q GW = 1.54 x 10" 14 , both 
models are safely below the nucleosythesis bound in Eq.(51). 

In the following we demonstrate the details of the resulting spectra h(k, rjn) and Q g (k) of RGWs, their explicit 
dependence upon the model parameters (3, (3 S , 7, and the modifications by the neutrino damping. 

Figure 4 gives the spectrum h(u,r]H) as a function of the frequency v without neutrino free streaming. To 
show the dependence upon the inflationary models, for the fixed r = 0.22, U\ = 0.75, (3 S = —0.3, we plot 
h(v,r]H) in three models of (3 = —1.8, —1.9, and —2.02. It is seen that h(v,r]H) is very sensitive to f3. A smaller 
[3 will generate less power of RGWs for all frequencies. The details of the spectrum is similar to that given in 
Refs.[3]. 

Figure 5 gives the comparison of the sensitivity curve of the ground-based interferometer LIGO with the spectra 
of = —1.8, —1.9, and —2.02 from Fig. 4. Here the vertical axis is the root mean square amplitude per root Hz, 
which equals to 

h(u) 

(52) 



Note that, relevant to LIGO, the frequency range is (10 ~ 10 4 ) Hz, which is not to be affected by the neutrino 
damping. Obviously from the plot, the LIGO I SRD [4] is yet not able to detect the signals of RGWs in the 
j3 = —1.8 model even with a very large ratio r = 2.2. Therefore, LIGO is unlikely to able detect the RGWs, as it 
currently stands. The Advanced LIGO with greatly enhanced sensitivity [4] will be able to put a tighter constraints 
on the parameters. 

Figure 6 is a comparison of the LISA sensitivity curve with the spectra from Fig. 4 in the frequency range 
(10 -7 , 10°) Hz. Although these frequencies are lower than that for LIGO, it is still not to be affected by the 
neutrino damping either. Assuming that LISA has one year observation time, which corresponds to frequency bin 
Au = 3 x 10~ 18 Hz (i.e., one cycle/year) around each frequency. Thus, to make a comparison with the sensitity 
curve, we need to rescale the spectrum h(u) in Eq.(44) into the root mean square spectrum h(u, Av) in the band 
Au [17], 

h(u,Au) = h(u)J — . (53) 
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Figure 5: Comparison of the spectra with the LIGO I SRD Goal sensitivity curve that has already been 
achieved by S5 of LIGO [4]. The vertical axis is the r.m.s amplitude per root Hz defined in Eq.(52). 




Figure 6: Comparison of the spectra with the LISA sensitivity curve [23]. The vertical axis is the r.m.s 
spectrum defined in Eq.(53). LISA will be able to detect the inflationary models of (5 = —1.8 and —1.9. 

LISA [23]. The plot shows that LISA by its present design will be able to easily detect the RGWs in the inflationary 
model of /? = —1.8. If the ratio r > 0.22, LISA will also be able to detect the inflationary models of f3 = —1.9. 
However, LISA is unlikely to be able to detect the model of j3 = —2.02. Here Fig. 5 and Fig. 6 also correct the 
mistake of Ref.[3], where improper comparison is made with the LIGO data and the LISA sensitive curve. As will 
be seen in the following, the neutrino free streaming practically affects only the spectrum in a frequency range of 
(10~ 16 ~ 10 -10 ), therefore, Fig. 5 and Fig. 6 are not to be changed by neutrinos. 

The influence of reheating stage on the spectrum is shown in Fig. 7. The spectra for three different values of 
Ps = 0.5, 0, —0.3 are given. It is clear that whereas the spectrum is almost unchanged by (3 S in the large portion 
of frequency range v < 10 7 Hz, a larger f3 s will damp the amplitude in a high frequency range 10 7 ~ 10 9 Hz. 
However, around v ~ 10 9 Hz the spectrum begins to increase considerably. This feature of RGWs in the GHz range 
is very interesting, as this high-frequency range of RGWs is the scientific gaol of some electromagnetic detecting 
systems, such as the one using a Gaussian laser beam [24], or a circulating microwave beam [25]. However, the 
predicted spectrum for the very high frequency range v > 10 10 Hz is not reliable, since the energy scale of the 
conventional inflationary models are less than 10 16 Gev. 




Figure 7: The reheating affects the spectrum h(v, rju) only in very high frequency range v > 10 7 Hz. 




Figure 8: The dependence of h(u,rjH) upon the dark energy Q\ in the accelerating universe. 

The influence of the dark energy on the spectrum h(v,rm) is demonstrated in Fig. 8, where Qa = 0.0, 0.7, and 
0.75 are taken respectively. Over the whole range of frequency 10~ 19 ~ 10 10 Hz, the amplitude of spectrum is 
altered by the presence of S7a, but the slope remains the same. In regards to the amplitude, firstly, the spectrum 
in a matter-dominant universe of Qa = is higher than those in an accelerating universe with Qa > 0, as Fig. 
8 shows, roughly by a factor ~ 1.3. This feature is due to the fact the scale factor a{tu) in the accelerating 
model is greater than that in the matter-dominant models, as has been explained at the ending paragraph of 
section 3. Secondly, among the accelerating models, by an analysis of the expression of h(u,r]H), the amplitude 
of h{y,r)n) is proportional to (f^y) 3 = 1/(1 + ze) 3 , as has been explicitly shown in Ref.[3]. This phenomenon 
occurs basically because, starting from the time r]E up to the present time r]n, the scale factor a{j]) increases by a 
different amount in models of different 17a ; thus, stretching of the physical wavelengths and damping of the mode 
hk{rf) are different correspondingly. In the accelerating models with p\ being constant, one has approximately 
( afagj ) 3 — ^m/^Ai thus, the amplitude of h(u, oc Q m /Q\, i.e., the model with more dark energy component 
has relatively a lower amplitude of h(v,r}i{). Note that, in interpreting this relation h{y,r\n) oc Q, m /QA, the 
dark energy Oa should be large enough, say, $7a > &m< t° ensure the sufficiently accelerating expansion. This 
phenomenon is verified now in Fig. 8, for example, the amplitude of the Oa = 0.75 model over that of the Oa = 0.7 
is found to be (|f|)/(§^) « 0.8. 
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Figure 9: The neutrino free-streaming reduces h(v,r]H) in the frequency range 10 17 ~ 10 10 Hz. 




Log V [Hz] 

Figure 10: The spectral energy density Sl g (v) for various inflation parameter f3. 

Presented in Fig. 9 is the modification of the spectrum h(v,r)jj) by the neutrino free-streaming up to the first 
order approximation to Eq.(34). The effect is pronounced in the low frequency range (10~ 16 ~ 10 -10 ) Hz, where 
the amplitude of h(u,r]H) is reduced by a factor ~ 20% in comparison with the model without neutrino free 
streaming. Our analytical spectrum qualitatively agrees with the numerical result in Ref.[16] in the relevant range. 
As we have mentioned earlier, LIGO and LISA operating around ~ 10 2 Hz and ~ 10 -3 Hz, respectively, will not 
be able detect this neutrino damping. But CMB anisotropics and polarization may be affected by that. 

The dependence of spectral energy density O g in Eq.(47) on the inflationary models is illustrated in Fig. 10. 
For the purpose of clarity, the neutrino free streaming is not taken into account. Clearly, Q g is very sensitive to the 
parameter (5, and a larger (5 gives a higher Q, g . The model of (3 = —1.8 has an Q g too high, and as mentioned in 
paragraph before Eq.(51), it has already been ruled out by the nucleosynthesis bound of Eq.(51). The Advanced 
LIGO [4] will be able to detect RGWs with Q, g h 2 > 10~ 9 at v ~ 100Hz, and it might impose stronger constraints 
on other inflationary models. 

The impact of dark energy Q,\ on the spectral energy density O g is plotted in Fig. 11 for the model = —2.02. 
It is clear seen that the accelerating expansion of the universe will cause a decrease of the amplitude of Q g over 
the whole range of frequencies, and a larger $7a gives a lower Q, g . Obviously, the effect due to the acceleration of 
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Figure 11: The spectral energy density Slg(v) for different A larger Q\ yields a lower $l g {y). 
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Figure 12: The neutrino free streaming reduces Q g {y) in the range 1CP 17 ~ lCP 10 Hz. 

expansion of the universe cannot be simply ignored. 

The clamping effect of neutrino free-streaming on the spectral energy density Q g (u) is illustrated in Fig. 12 for 
the model j3 = —2.02. The effect is mostly within the frequency range 10 -16 ~ 10~ 10 Hz, where the amplitude of 
0, g (y) drops visibly by a factor of ~ 36%. Correspondingly, the energy density of RGWs in Eq.(48) is now reduces 
to Slew = 1.1 x 10~ after considering the neutrino free-streaming. Recall that it was Qgw = 1-54 x 1CP 14 
without neutrino damping, thus the neutrino damping has caused a drop of ~ 29% of £Igw- Our result is also 
qualitatively consistent with the numerical calculation in Ref. [16]. 
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Appendix 

In this appendix we derive the anisotropic stress tensor Kij of cosmic neutrinos during their free streaming, 
filling in the details skipped in Ref.[14]. Next we present the perturbation method to systematically solve the 
equation of GRWs. As a merit, this method applies for the realistic situation of a time-dependent fractional 
enercv Hensitv f..(n\ nf neutrinos a nnti7Prn Herminlincr time 7/.j_„ zk fl and a nnn-vaniqhincr time derivative 



x'( u dec) 7^ 0. This solution is an extension to the previous works in Refs.[14] [15]. The same notations as in 
Ref.[14] is used. 

In the radiation stage at temperatures > 2 Mev , the neutrinos are in equilibrium with other relativistic species, 
such as electron and photons, together forming the radiation component. During this period, practically all the 
fc-modes hk(rj) of cosmological interest are still far outside the horizon, thus remain a constant, hk(r)) ~ constant, 
thus, are not effected by the neutrinos. With the further expansion of the universe when the temperature drops 
down to < 2 Mev, the neutrinos are going out of equilibrium and are starting to stream freely in space. Then the 
neutrinos will be able to influence the modes hk(rj) of short wavelengths that re-enter the horizon. The equation 
of RGWs then becomes inhomogeneous 



h'l 3 { V ) + 2^^(r?) - V 2 hiM = l^Ga 2 ^M- (54) 



Here the source term 7Ty, contributed by neutrinos, is the anisotropic part of the stress tensor and is effective 
only during the period rj dec < V < V2< from the neutrinos decoupling up to the beginning of the matter domination. 
When the matter domination begins, the neutrino number density has been diluted out by a factor ~ 10~ 3x6 , so 
the source ir^ is effectively switched off after the matter domination. In terms of the neutrino distribution function 
n(x, p,t) and the momentum p l , the spatial part of the neutrino energy-momentum stress tensor is written as 

T) = -±=J d 3 pn(x, p, t)f Pj /p°. (55) 

To keep the same notation with Ref.[14], here the cosmic time t = J a(rj)dn is used. In the presence of the 
perturbations of the metric, n(x, p, t), p l , and p° all depends on hij. So the stress tensor is written as a sum 
of 

T j =P6 ij + Tr ij , i,j = 1,2,3 (56) 

where PSij is the unperturbed part with P being the homogeneous and isotropic pressure, and the anisotropic 
stress tensor ir^ is the perturbed part caused by hij 

7T^-(x,t) = -L<5 / d 3 pn(x,p,t)^. (57) 

Since is small, only the first order of hij is needed in evaluating tv^. The distribution function n(x, p, t) 
satisfies the Boltzmann equation 

dn dn dx l dn dpi 
~dt + d^lF + dp~^t =0 ' 
With dx % /dt = p l /p° and the geodesic equation dpi/dt = ^g^iP^p 6 = \hjk^p h , Eq.(58) can be expanded as 

dn p i dn 1 pip k dn 

m + ^ + 2 h ^^rdp- = - m 

At the instant t^ ec of the decoupling the neutrinos are still an ideal gas, so one writes 

n(x, p, t) = n(x, p, t dec ) + <5n(x, p, t), (60) 

where 



N r / / \ 1-1 

n(x, p, t dec ) = ——3 [ exp [y/gv ( x , t dec )piPj/T dec j + lj (61) 



(2vr) 2 



is the distribution function of the ideal gas of temperature r^ ec , and Sn represents the perturbation satisfying 
Sn = at tdec- In our treatment pi is treated as the unperturbed momentum and p l = g^pj as the perturbed 
one. Substituting Eq.(60) into Eq.(59), neglecting the higher order term ^jfc,* 2 ^-^ 1 . using dnd ec /dt = 0, 

= -\n{p) ppjPk-^-hjk^, tjec), (62) 

(Ref.[14] missed a factor \ in Eq.(62)), and dn/dpi = n'{p)pi, where n'{p) = dn/dp, pi = pi/p and p = ^PiPi , 
the Boltzmann equation reduces to the following 

dSfi p' dSfi p d 

~~dT + W)^ = -Wf' mp ^ k ^ {hij{ ^ " M *' (63) 

where n(p) = j^j[ex.p(p/Tde C adec) + I] -1 - This equation can be decomposed into the Fourier k-modes, as in 
Eq.(19), and each mode has the formal solution 

Sn k (p, 7]) = - l -pn(p)piPjp • k [ dij'jP^'-ri \hij(k, rf) - h i:j (k, %ec )] , (64) 

where the comoving time dw = dt/a is used, and the integrand function depends on hij explicitly, where 

h tJ (k, V )^Y,^ h k\v)- (65) 

a 

Defining the variable u = k(r] — % ec ), Eq.(64) just reduces to Eq.(13) in Ref.[14]. There are four terms that 
contribute to ttij in Eq.(57). Specifically, the perturbations to the distribution function n give two terms: 



/<fofc(% ec ) = -7T n ' (P)PiPjPmPlh m l(k,r} dec ), (66) 



p0 2a 

n2 



^Sn k ( V ) = -^n'WpiPjtinPl \hml(k, v) ~ h m i(k, r, dec ) - f drj'e^'^hUk, v')] , (67) 
where integration by parts with respect to the time rf has been used, the following two terms also contribute 

flip) ■ J) 

—— Pj Sp t = --fipip j h u (k,r]), (68) 

~'^M Sp ° = ^ fi PiPjPmPih m i(k,rj). (69) 

One puts these four terms from Eq.(66) through Eq.(69) into Eq.(57) and carries out the integration J d 3 p. The 
spherical coordinates with z = k can be used in doing the angular integration, so that p ■ k = cos# = fi. Since 
h m i is transverse, one has 

J dVLpiPjp" m pih mi = jiSimSji + 5u5 jm )h m i J ^ - fi 2 ) 2 . (70) 

After some calculation, one arrives at the resulting anisotropic stress 

ni = -4P»(V) P dn'K{kn - krf^irf), (71) 

where p v (rj) = a' 4 j d 3 ppn(p) is the neutrino density, and K is the kernel defined as 

K(s) = 1 £ d^ (1 - M 2 ) 2 e^ = (72) 



with j2(s) is the spherical Bessel function. Since hij is traceless and transverse, so is ir^, by Eq.(71), with n l i = 
and nijj = 0. Substituting Eq.(71) into Eq.(54) and using the Friedmann equation a' 2 /a A = 8irGp/3 yields the 
integro-differential equation [14] 

h'Liri) + ^Mh' k ( V ) + k 2 h k (rj) = -2Af v {i 1 )\ , ^h 2 f drj'K(k V - krf)h' k (rf), (73) 

a vl) a {V) Jri dec 

where the fractional neutrino energy density 

/.(,) = $$■ (74, 

Although at present the dark energy Q\ is dominant, but it is negligible during the radiation-dominant stage, in 
comparison with the matter, neutrino, and radiation components. Even in the dynamic models of dark energy 
evolving with time, the contribution from pa(v) during radiation-dominant stage is not allowed to be more than 
a few percent of the total energy [26] [27]. Therefore, Eq.(74) is practically equal to 

fAn) = ; Jfl , (75) 
1 + a(ri)/a e q 

where a eq = 0(772) is the scale factor at the radiation-matter equality, and 

'"<°> = STS- (76) 

with Q u and J7 7 being the present fractional energy density o f the neutrinos and the radiation, respectively. Since 
Vdec ^ ?7e. we can write a(r]) = a e r] with a e defined in Eq.(4) for the radiation-dominant stage. Introducing the 
variable u = krj and setting 

hk(u) = h k (r) dec )x(u), (77) 

then Eq.(73) is reduced to 

X » + - X '(u) + X {u) = -24 fA f ] f dUK(u - U) X '(U), (78) 
u u z {l + au) J Udec 

where a = a e /ka eq . 

Dicus ii Repko [15] present an analytical solution of Eq.(78) under the approximation of setting Ud ec = 0, 
a = and x'iVdec) = 0. This is only valid for those modes that reenter the horizon well after the neutrino 
decoupling. Note that, in the coordinate that is used, the decoupling time Ud ec / as in Eq.(10). Besides, for 
the short wavelength modes the derivative x'iVdec) 7^ 0. Moreover, actually au = 1 at the time 772, thus, setting 
a = for the whole period (7/^,772) would lead to an over-account of the fractional neutrino energy density 
fu(v) ar| d consequently would give a lower amplitude of RGWs. Unlike in Refs.[14, 15] , here we do not make the 
above-mentioned approximation, but instead, keep Ud ec < a and x'iVdec) as they are. We use a perturbation method 
to solve the integro-differential equation (78) analytically, which can achieve high accuracy as one requires. Note 
that the source term on the r.h.s. of Eq.(78) is relatively small, and setting it to be zero yields the homogeneous 
equation 

XoH + ^XoH+XoH=0, (79) 

with the solution 

e iu e~ iu 

Xo (u) = c l — + c 2 , (80) 



as the th order approximation to Eq.(78), where c\ and c 2 are the coefficients determined by the continuity 
condition at Ud ec - One substitutes x'o( u ) i n place of x'( u ) i n the integration of Eq.(78) to give the I s * order 
approximation 

X'l(u) + -x'M + xi in) = -24 f| 0) f dt^u - U) X ' (U), (81) 
u u 2 (l + au) J Udec 

which is a differential equation with a known inhomogeneous term. It has a particular solution 

X{) L c W[ yi ,y 2 ](v) r{V)dV 

24/, (0) r sin(n - v) r> j 2 (v - s) , 

= ~— L dv WT^) L ds ¥^¥ Xo{s) > (82) 

where r(v) represents the inhomogeneous term of Eq.(81), and yi = ^j- , U2 = are the two linearly independent 
solutions to the homogeneous counterpart, and W[yi,y2]{v) = —2i/v 2 is the Wronskian. Therefore, the solution 
of Eq.(81) is given by 

Xi(u)=xo(n)+X*(u), (83) 

which is also the I s * order approximate solution of Eq.(78). Similarly, one substitutes the I s * order solution xi{ u ) 
into Eq.(82), and obtains the 2 nd order approximate equation, 

X ' 2 '(u) + -x' 2 {u) + X2{u) = -24 J;f ] [ U dUK(u - U) X [(U), (84) 
u u\l + au) J Udec 

which has a particular solution 

**/ \ 24/,(r? dec ) f u sm(u-v) f v j 2 {v - s) , 

X («) = / dv 77~~ r / ds- ro-Xi(s), 85 

u Ju dec v(l + av) Ju dec (v - s) 2 

and thus the 2 nd order approximate solution is 

X2(u)=Xo(u)+X*»- (86) 

By the same routine, the higher order solutions can be obtained. In fact, as our calculation shows that the 1 st order 
approximation is already accurate enough for the purpose of computing the spectrum for RGWs. An important 
advantage of our solution is that it is valid for those modes that reenter the horizon before or after the decoupling 
time rjdec- Integrations, such as in Eqs.(82) and (85), can be done easily by common computing tools. 
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